% Folder containing the color-coded labeled cluster images
folderPath = '';  % <-- CHANGE this to your folder

% Image file type (change to *.jpg if needed)
imageFiles = dir(fullfile(folderPath, '*.png'));

% Output folder
resultsFolder = fullfile(folderPath, 'Filtered_By_Circularity');
if ~exist(resultsFolder, 'dir')
    mkdir(resultsFolder);
end

% Excel file path
excelFilePath = fullfile(resultsFolder, 'Filtered_Cluster_Data.xlsx');

% Filtering thresholds
minClusterSize = ;
maxCircularity = ;

for i = 1:length(imageFiles)
    try
        % Load image
        imageName = imageFiles(i).name;
        imagePath = fullfile(folderPath, imageName);
        rgbImage = imread(imagePath);

        % Reconstruct binary mask from RGB image (assumes white background)
        grayImage = rgb2gray(rgbImage);
        binaryImage = grayImage < 255;  % clusters are colored, background is white

        % Initial labeling
        labeledImage = bwlabel(binaryImage);
        stats = regionprops(labeledImage, 'Area', 'Perimeter', 'Centroid', 'PixelIdxList');

        % Compute circularity
        circularity = zeros(length(stats), 1);
        for j = 1:length(stats)
            if stats(j).Perimeter > 0
                circularity(j) = 4 * pi * stats(j).Area / (stats(j).Perimeter^2);
            end
        end

        % Filter based on area and circularity
        validClusters = [];
        for j = 1:length(stats)
            if stats(j).Area >= minClusterSize && circularity(j) < maxCircularity
                validClusters = [validClusters; j];
            else
                labeledImage(stats(j).PixelIdxList) = 0;  % remove invalid
            end
        end

        % Rebuild cleaned binary image from valid clusters
        filteredBinary = labeledImage > 0;

        % Relabel for display and measurement
        labeledClean = bwlabel(filteredBinary);
        statsFiltered = regionprops(labeledClean, 'Area', 'Centroid', 'Perimeter');

        % Skip image if no clusters left
        if isempty(statsFiltered)
            disp(['No valid clusters in image: ', imageName]);
            continue;
        end

        % Get shape properties
        areas = [statsFiltered.Area]';
        perimeters = [statsFiltered.Perimeter]';
        centroids = cat(1, statsFiltered.Centroid);
        circularityFiltered = 4 * pi * areas ./ (perimeters.^2);
        clusterIDs = (1:length(areas))';

        % Generate display image
        numClusters = length(clusterIDs);
        uniqueColors = distinguishable_colors(numClusters, [1 1 1]);
        rgbOut = label2rgb(labeledClean, uniqueColors, 'w', 'shuffle');

        % Overlay IDs
        figure; imshow(rgbOut); hold on;
        for k = 1:numClusters
            text(centroids(k,1), centroids(k,2), num2str(clusterIDs(k)), ...
                'Color', 'black', 'FontSize', 14, 'FontWeight', 'bold', ...
                'HorizontalAlignment', 'center');
        end
        hold off;

        % Save annotated image
        [~, fileName, ~] = fileparts(imageName);
        outputImagePath = fullfile(resultsFolder, [fileName '_filtered_by_circularity.png']);
        saveas(gcf, outputImagePath);
        close;

        % Save table to Excel
        clusterData = table(clusterIDs, areas, centroids(:,1), centroids(:,2), perimeters, circularityFiltered, ...
            'VariableNames', {'ClusterID', 'Area', 'Centroid_X', 'Centroid_Y', 'Perimeter', 'Circularity'});

        sheetName = fileName(1:min(31, length(fileName)));
        writetable(clusterData, excelFilePath, 'Sheet', sheetName, 'UseExcel', false);

        disp(['Processed and saved: ', imageName]);

    catch ME
        warning(['Error processing ', imageName, ': ', ME.message]);
    end
end

disp('All images have been processed and saved based on circularity filtering.');

% --- Color function ---
function colors = distinguishable_colors(n_colors, bg)
    if nargin < 2
        bg = [1 1 1];
    end
    colors = lines(n_colors);
    if n_colors > size(colors, 1)
        colors = rand(n_colors, 3);
    end
end
